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Abstract We demonstrate the remarkable effectiveness of boundary value formulations coupled to numerical contin- 
uation for the computation of stable and unstable manifolds in systems of ordinary differential equations. Specifically, 
we consider the Circular Restricted Three-Body Problem (CR3BP), which models the motion of a satellite in an Earth- 
Moon-like system. The CR3BP has many well-known families of periodic orbits, such as the planar Lyapunov orbits and 
the non-planar Vertical and Halo orbits. We compute the unstable manifolds of selected Vertical and Halo orbits, which 
in several cases leads to the detection of heteroclinic connections from such a periodic orbit to invariant tori. Subsequent 
continuation of these connecting orbits with a suitable end point condition and allowing the energy level to vary leads 
to the further detection of apparent homoclinic connections from the base periodic orbit to itself, or the detection of 
heteroclinic connections from the base periodic orbit to other periodic orbits. Some of these connecting orbits are of 
potential interest in space mission design. 

Keywords Restricted three body problem • Boundary value problems • Invariant manifolds • Connecting orbits • 
Numerical continuation 



1 Introduction 

Numerical continuation of solutions to boundary value problems (BVPs) has been used extensively to study periodic 
orbits and their bifurcation s, including homoclinic and heteroclinic orbits, in a wide variety of systems. For a recent 
overview see the articles in lKrauskopf et al.l J2007I) . In particular, these techniques have been applied to compute f am- 
ilies of periodic orbits in the Circular Restricted 3-Body Problem (CR3BP); see, for example, iDoedel et al.l (120031) . In 
recent years inva riant manifolds such as those in the Lorenz system have also been c omputed in detail using numeri- 
cal continuation ( Krauskopf and Osinga 1 120071 : lAguirre et alJEoTTI : IDoedel et alJEoTTh . At the same time, continuation 
methods have been developed for computing and continuing homoclini c and heteroclinic connectin g orbits between pe- 
riodic orbits and equilibria or periodic orbits (IDoedel et all 2008, 20091; iKrauskopf and R ieB 2008). In the current work 
we use a combination of these techniques to illustrate their effectiveness in computing stable and unstable manifolds 
and connecting orbits in the CR3BP. 

Th e re is much l i teratu r e on invariant manifo l ds and connecting orbits i n the CR3BP; see for ex ample iGomez et al.1 
d2QQ4h : iKoon et alJ (120081) ; ILo and Rossi dl997l) ; iDavis et al.l ( bold . 1201 lb : iTantardini et all d201Ch . In particular, con- 
necting orbits in the planar CR3BP are well un derstood. The existence of connecting orbits in the planar probl e m has 
been proved analytically in lLUbreetalJdl985h . and by computer assisted methods in Iwilczak and Zgliczvhskil (I2003L 
l2005h . Furthermore, these orbits have bee n extensively stu died numerically using initial-value techniques and semi- 
analytical tools; see iBarrabes et al.l d2009r) : ICanaliasI d2006r) and references therein. In the case of initial-value tech- 
niques the initial conditions are varied in order for an appropriately chosen end point condition to be satisfied. This 
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approach is commonly referred to as a "shooting method" and, for a more stable version, "multiple shooting". Initial- 
value techniques can also be very effective in the computation of invariant manifolds in the CR3BP. However, sensi- 
tive dependence on initial conditions may leave parts of the manifolds unexplored, unless very high accuracy is used . 
Other efficient methods for computing invaria nt manifolds include semi-analytical approximations (Ijorba et al .11 19991 : 
lAlessi et aDl2QQ9l : iGomez and Mondeloll200ll) . The latter methods are very precise in a ne ighborhood of the c enter of 
expansion, and rely on other methods to extend the manifolds outside these neighborhoods jGomez et al.ll200ll). 



Invariant manifold techniques around libration points have been used successfully in mission design ( Lo et al. 2004) 



The Genesis spacecraft mission, designed to collect samples of solar wind and return them to the Earth (ILo et al .112001 



|Lo et al.| 2 
(ILo et al.l l2 



is often considered as the first mission to use invarian t manifolds for its planning, while other missions have used 
libration point techniques dDunham and Farquharll2003l) . Having a precise idea of the geometry of invariant manifolds 
and their connections is desirable in the design of complex low thrust missions. 

Using our continuation approach we construct a continuous solution family in the manifold as the initial value is 
allowed to vary along a given curve. The continuation step size governs the distance between any computed trajectory 
and the next trajectory to be computed. Here "distance" includes the change in the entire trajectory, and not only in 
the initial conditions. In fact, this distance typically also includes other variables in the continuation process, such as 
the integration time and the arclength of the trajectory. This formulation allows the entire manifold to be covered, up 
to a prescribed length, integration time, or other termination criterion, even in very sensitive cases. Special orbits, such 
as connecting orbits to saddle-type objects, can be detected during the continuation. For example, a straightforward 
computation of the Lo renz manifold in thi s fashion yielded up to 512 connecting orbits having extremely close initial 
values, as rep orted i n Doedel et al. | (|200d^ In related work, the intersections of the Lorenz manifold with a sphere 
are studied in lAguirre et al.l ( 201 lh and Doedel et In the case of the Lorenz manifold, the sensitivity on 

initial conditions results from the significant difference in magnitude of the two real negative eigenvalues of the zero 
equilibrium that give rise to this manifold. Fixed precision integration in negative time may cover only part of this 
manifold, missing a portion in and near the direction of the eigenvector of the smaller negative eigenvalue. 

In this paper we give an overview of continuation techniques, as used to compute periodic orbits, invariant manifolds, 
and connecting orbits. We also give several examples that illustrate how these techniques provide an effective and 
relatively easy-to-use tool for exploring selected portions of phase space. The richness of the solution structure of the 
CR3BP limits the extent of our illustrations. However, the techniques presented her e are expected to be useful in further 
studies. In this respect, the current version of the freely available AUTO software tooedel et alJbOlCh includes demos 
that can be used to re-compute some of the numerical results presented here, including their graphical representation. 
These demos can also be adapted relatively easily to perform similar numerical studies of stable and unstable manifolds 
of other periodic orbits in the CR3BP that are not considered in this paper, as well as for entirely different applications. 

This paper is organized as follows. In Sect. [2 we recall some well-known facts about the CR3BP, namely, its equi- 
libria, the libration points, and the basic periodic solution families that will be considered in this paper, namely the 
planar Lyapunov orbits and the Vertical, Halo and Axial families. In Sect. [3] we review how boundary value techniques 
are used to compute periodic orbits in conservative systems, and how these techniques can also be used to compute the 
eigenfunctions associated with selected Floquet multipliers. These data provide a linear approximation of the unstable 
manifold of the periodic orbit. 

Sect. HI describes the continuation method used for computing unstable manifolds of periodic orbits. This involves 
first setting up an extended system with both the periodic orbit and its eigenf unction. Using the resulting information 
an initial orbit within the manifold is computed. This orbit is then continued, as its starting point is free to vary along a 
line that is tangent to a linear approximation of the unstable manifold, thereby tracing out the manifold. The algorithm, 
using pseudo-arclength continuation, is not guaranteed to compute the whole manifold in a single computation, because 
obstacles may be encountered during the continuation. In such cases the manifold may be completed, for example, by 
ad ditional continuation s from different starting orbits, or by using a suitably adapted continuation procedure as is done 
m iDoedel et alJ feOllh . However, the obstacles themselves are also of interest. They may correspond to orbits in the 
unstable manifold that require an arbitrary long time interval to reach a specified termination plane, because they pass 
arbitrarily close to a connecting orbit between the original unstable periodic orbit and another invariant object. In this 
way c onnecting orbits can be detected, as was done to detect the 512 heteroclinic connections presented in lPoedel et al.1 
J2006h . 

In Sect. [5] we show the results of computations of unstable manifolds of Vertical Vi and Halo Hi orbits. When the 
manifolds are computed to a sufficient distance from the original periodic orbit, we find what appear to be heteroclinic 
connecting orbits from the original periodic orbit to an invariant torus. The tori found this way must have saddle-type 
instability, since the connecting orbit approaches it, but ultimately also leaves the neighborhood of the torus. Such 
connecting orbits may be more difficult to find with initial value integration. 

In Sect. [6l we describe a method for continuing the periodic-orbit-to-torus connections as solutions of an 18- 
dimensional ODE, when the energy is allowed to vary. These computations lead to the detection of other interesting 
connecting orbits. Sect. [7] discusses three representative examples. First we consider a family of connections from Hi 
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Halo orbits. These connections loop once around the Earth before approaching an invariant torus, which is itself close to 
the Hi orbit. We refer to such a torus as a quasi-Hi torus. During the continuation, with changing energy of the originat- 
ing Hi Halo orbit, we encounter a number of interesting connecting orbits. Specifically, we find a homoclinic orbit from 
an Hi Halo orbit to itself that loops once around the Earth, a heteroclinic connection from a northern Hi Halo orbit to its 
southern counterpart, a heteroclinic connection to a planar Li Lyapunov orbit, and a connection to a 5:1 resonant orbit 
on a torus near the corresponding southern Halo orbit. We also find a connecting orbit from an Hi Halo orbit to a torus 
on which the orbit bounces back and forth between a northern Axial Ai orbit and its southern counterpart. Each of these 
special connecting orbits occurs for a specific energy of the originating Hi Halo orbit (and of course, by conservation 
of energy, the orbit it connects to has the same energy). Secondly we study a family of connecting orbits from an Hi 
Halo orbit that loop four times around the Earth. We find connections to an L2 Lyapunov orbit, an H2 Halo orbit, and 
to a 5:1 and a 6:1 resonant torus near the libration point J^2- Thirdly, we consider a family of connecting orbits on the 
Moon-side of the unstable manifold of the Hi Halo orbits that connect directly (without looping around the Earth) to a 
torus near «£?2. We find an example of a direct connecting orbit from an Hi Halo orbit to a planar L2 Lyapunov orbit. 
To the best of our knowledge, these connecting orbits have not been found before. We must stress, however, that from a 
space mission design point of view, these orbits are sensitive to initial conditions and would require control techniques 
to stay on them. 

In Sect. [U we discuss global theoretical aspects of our results and their relation to the existing literature. In partic- 
ular we see that the connecting orbits from Hi Halo orbits to Ai Axial orbits or to Li or L2 planar Lyapunov orbits 
are codimension-one in the dynamical systems sense, and hence should occur for specific values of the energy of the 
originating Hi Halo orbit, as we observed numerically. In contrast, homoclinic and heteroclinic connecting orbits be- 
tween Hi Halo orbits, which were observed numerically, are codimension-two, and so should not normally occur. We 
show that the connecting orbits from the Hi Halo orbits to quasi-Hi tori are generic, and suggest that the numerically 
observed homoclinic and heteroclinic connecting orbits between Hi Halo orbits are actually connections to quasi-Hi 
tori where the minor radius of the torus is so small as to make the torus visibly (and for the purpose of space mission 
design) indistinguishable from the Hi Halo orbit that it envelopes. 

Finally in Sect. [9] we discuss some computational aspects of our numerical computations, such as the discretization 
used, and typical computer time needed. 



2 The Circular Restricted 3-Body Problem 



The CR3BP describes the motion of a satellite S with negligible mass in three-dimensional physical space. The motion 
is governed by the gravitational attraction of two heavy bodies, which are assumed to rotate in circles around their 
common center of mass; see Fig.[TJa). In this paper we call the heaviest body the "Earth" E, and the other heavy body 
the "Moon" M, and we use their actual mass ratio, namely, /1 w 0.01215. Any other mass ratio is allowed, such as for 
the Sun- Jupiter system, with a mass ratio of ji w 0.0009537. Without loss of generality the total mass can be scaled to 
1, so that the Earth and Moon have mass 0.98785 and 0.01215, respectively. 
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The equations of motion of the CR3BP as given in lDanbvl (1 19921) are 

x = 2y + x - (1 -fi) — 3 (I 3 — 



y= -2i + v-(l- M )^-M^ 
where is the position of the zero-mass body, and where 



(1) 



(x + n) 2 +y 2 +z 2 , r 2 = y(jc- 1 + ju) 2 + y 2 +z 2 , 

denote the distance from S to the Earth and to the Moon, respectively. The CR3BP has one integral of motion, namely, 
the energy E: 



x 2 +f+i 2 



+ U(x,y,z) , U(x,y,z) = ~\(x 2 +y 2 ) - l -JL _ E _ j 



where C/ (jc, v,z) is the effective potential. Astronomers also often use the Jacobi constant C, defined as C — —2E. 
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Fig. 2: Bifurcation diagrams of families of periodic solutions of the Earth-Moon system bifurcating fro m the libration 
points (a) and (b) J^2- A detailed description of the families represented in this diagram can be found in lDoedel et al.l 
d2QQ7h . The "thick" portions of the curves labeled Hi and H2 denote periodic orbits where all 6 Floquet multipliers 
are on the unit circle: 1, l,exp(±/c),exp(±/d). For thin solid portions exactly two real Floquet multipliers are off the 
unit circle: 1 , 1 , a, l/a,exp(±«c). The dotted portions denote periodic orbits where all 6 Floquet multipliers are real: 
1, l,a, l/b. Here a,b,c,d G R, \a\ > 1, \b\ > 1, and c,d E (0, 7r). The small squares labeled Ly and Vn denote 

branch points. 



Libration points, in both the planar and three-dimensional spatial system, are equilibrium points in a co-rotating 
frame; see Fig. [Tib). As already used, we denote the libration points by Jzfi, • • • , J^. There are families of periodic 
orbits (in the co-rotating frame) that bifurcate from each of these libration points, and we refer to these as the primary 
families. Many more families subsequently bifurcate from the primary families. We refer to the bifurcation points as 
branch points. Several fa milies of periodic solutions of the Earth-Moon system are represented in Fig. |2j see also 
iDoedel et al.l d2003ll2007h and references therein. In the present work we focus on four families, namely, the Vertical 
orbits Vi, the planar Lyapunov orbits Li, the Halo orbits Hi, and the Axial orbits Ai. The families of planar Lyapunov 
orbits Li and the Vertical orbits Vi emanate directly from the libration points JzfJ; these are primary families (Fig. [3}. 
The families of Halo orbits bifurcate from the families of Lyapunov orbits Li, while the families of Axial orbits connect 
and bifurcate from the Vi and Li families. We refer to the Halo and Axial orbits as secondary families (Fig. |4j). These 
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families are all well-documented in the literature, but their names are som etimes di f ferent . For example , the Halo, Axial, 
and V ertical orbits are known as type "A", "B", and "C", respectively, in lGoudas" I dl96lh and lHenon I (11973k Farquhar 

1968 ) coined the name "Halo" for that family. The term 'Axial" comes from Dpedel et al] J2007I) , whereas iDoedel et al 

2003) used the term "Y" for "Yellow". 




Fig. 3: (a): A selection of periodic orbits from the planar Lyapunov family L\, which bifurcates from the libration point 
Jzfi, as seen in the bifurcation diagram in Fig.[2ja). Labeled are the special Lyapunov orbits Ln from which the Halo 
family Hi emanates, and L12 from which the Axial family Ai bifurcates. In this and subsequent figures the small cubes 
denote libration points, (b): Selected orbits from the Vertical family Vi, which also bifurcates from the libration point 
S£\ . In the linear approximation near ££\ these orbits are indeed "vertical", i.e., x and y are constant along it, with y = 0, 
while z oscillates around zero. The Axial family bifurcates from the orbit labeled Vn. 



3 BVPs for Periodic Orbits and Eigenfunctions 

In this and the next section, we describe the algorithms used in the various stages of the computations in some detail, 
so that it will be possible for the reader to replic ate the algorithms i n other applications, and to better understand the 
functioning of the downloadable CR3BP demos (IDoedel et al .1120 101) . Specifically, in this section we explain the pre- 
liminary computations that precede the actual computation of stable and unstable manifolds of periodic orbits, namely, 
the computat i on of the pe riodic orbits themselves and of their associated eigenfunctions. The discussion follows that in 
IDoedel et alJ fe00il2008h . 

To formulate a suitable boundary value problem in AUTO, the second order system of ODEs (0]) first needs to be 
rewritten as a six-dimensional first order system, 

u(0 =i(u(t),n), f:M 6 xR^R 6 , 

where u = (x, y,z, v x , v y , v z ) = (x,y,z,x,y,z). As usual, when we plot the orbits graphically we project into M? by only 
plotting the spatial coordinates (x,y,z). Time is scaled to the interval [0, 1] in the BVP formulation for computing a 
periodic orbit, which changes the system of differential equations to 

u(t) = rf(u(0,M), 

where T is the period of the periodic orbit. In addition, for a conservative system with one conserved quantity, we need 
to add a term with an "unfolding parameter" in order for the BVP continuation computations to be formally well-posed 
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(b) 



Fig. 4: (a): Selected orbits from the Halo family Hi that bifurcates from the Lyapunov family Li at Ln. (b): Selected 
orbits from the Axial family Ai that connects to the Lyapunov family Li at L12 and to the Vertical family Vi at Vn : see 
also Fig. 12 



JPoedel et al.ll2007l ; iMunoz-Almaraz etaD 12003): see also iMunoz-Almaraz et al.l (120071) . A suitable and convenient 
choice in the specific case of the CR3BP is the term crd(u), where d(u) = (0,0,0, v x ,v y ,v z ). The vector field with the 
unfolding term then becomes 

u(0 = rf(u(r),M) + ^d(u(0), 

which from here on we simply write as 

u(0 = rf(u(0,a), (2) 

also omitting the mass-ratio parameter /1, as it is typically fixed in the computation of families of periodic orbits. Notice 
that the specific choice of unfolding term d(u) used here would represent a damping (or forcing) term if o ^ 0, which 
would preclude the existence of periodic orbits. However, the unfolding parameter o is one of the unknowns in the 
continuation procedure, and will always be zero (to numerical precision) once solved for. Thus the unfolding term is 
simply a technical device necessary to obtain well-posedness of the BVP, and we do not force or damp the equations of 
motion. 

Written in full, the system is therefore given by 

x = T v x , 
y = T v y , 
z = T v z , 

o a (3) 

v x = T [2v ); +x-(l- i u)(x + i u)r 1 3 - i u(x-l + i u)r 2 3 ] + av x , 
v y = T [-2v x +y-(l-Li)yr- 3 -!Ayr- 3 } + ov y , 
v z = T[-(l- n)zrr 3 - \izr 2 - 3 \ + gv z . 
To complete the BVP formulation we need to add the periodicity equations 

u(l) = u(0) . (4) 

If u(t) solves Eqs. ^ and © then u a (t ) — u(t -j- oc) i s also a solution for any time-translation a. To specify a unique 
solution we impose the phase constraint fooedellll98l1) 



1 



\u{t)Mt))dt = 0, (5) 
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where uo(f) is the preceding computed solution along the solution family. Furthermore, for the purpose of continuing 
a family of periodic solutions, we add Keller's pseudo-arclength constraint (Keller ! 1977b . which in the current setting 
takes the form 

/ (u(0-uo(0»wo(0>* + (r-2b)rj + (a-a )o^ = As, (6) 
Jo 

where (uo, To, Ob) corresponds to a computed solution along a solution family, (u, T, a) is the next solution to be com- 
puted, and As is the continuation step size. The notation " ' " denotes the derivative with respect to As at As — and (•, •) 
denotes the dot product. Since we are dealing with a conservative system, we already have families of periodic solutions 
even when the mass-ratio /I is fixed. For a computed solution (uo, To, Ob), and a given step As, the unknowns to be 
solved for in any continuation step are the periodic orbit u(t), its period T, and the unfolding parameter a. Eq.[6]forces 
all of these unknowns to be close to those of the previous solution. In particular, u(f ) must be close to uo(f) for all t, and 
not just for t = 0. We reiterate that the unfolding parameter o is an active unknown in the computations that regularizes 
the boundary value formulation, although it will be found to be zero to numerical precision. During the continuation 
the Floquet multipliers of the periodic solution are monitored by computing a special decomposit ion of the monodromy 
matri x that arises as a by-product of the decomposition of the Jacobian of the collocation system dFairgrieve and JepsonI 
E99l|). 

For the purpose of computing a stable or unstable manifold of a periodic orbit we need the corresponding Floquet 
eigenf unction. Specifically, we assume that the periodic orbit has a single, real, positive Floquet multiplier outside the 
unit circle in the complex plane, which gives rise to a two-dimensional unstable manifold of the periodic orbit in phase 
space. The eigen function correspond ing to this multiplier provides a linear approximation to the manifold close to the 
periodic orbit. In lDoedel et al.l (120081) it is shown that this eigenfunction can be obtained as a solution v(f ) of the BVP 

v(0 = rf u (u(0,o)v(0+Av(0, 

v(l) = ±v(0), (7) 
(v(0),v(0))=p, 

where v(l) = +v(0) in the case of a positive multiplier, and v(l) = — v(0) in the case of a negative multiplier. Here, A 
is the characteristic exponent and the corresponding Floquet multiplier is given by ±e^ . Eq. (0 represents the Floquet 
eigenfunction/eigenvalue relation for the linearization of Eq. © about a periodic orbit u(f). The norm of the value of 
the eigenfunction at time t — is normalized to be y/p , where typically we use p = 1 . If only one Floquet multiplier 
is real and greater than one in absolute value, this gives a unique (up to sign) unstable eigenfunction \(t). Likewise, a 
unique stable eigenfunction is obtained if only one Floquet multiplier is real and less than one in absolute value. In our 
illustrations we only compute unstable manifolds, so we have A > in Eq. (0; however all algorithms apply equally 
well to stable manifolds. We also restrict to the case where the Floquet multiplier of interest is positive, so v(l) = +v(0) 
in Eq. ©, and the corresponding manifold is orientable rather than twisted. A linear approximation of the unstable 
manifold at time zero is then given by 

r(0) = u(0) + ev(0) , (8) 

for £ small. 

An alternative formulation is to put the actual Floquet multiplier in the boundary condition rather than in the lin- 
earized differential equation, using the variational equation \(t) = rf u (u(f),0)v(f) and boundary condition v(l) = 
fcv(0), where K is the actual Floquet multip lier. However, the f ormulation in Eq. ©, as used here, has been found to 
be more appropriate for numerical purposes tooedel et alJ l2008\ This is related to the fact that the multipliers, i.e., the 
values of K — e^, can be very large or very small. 

The algorithmic steps that lead to the linear approximation of the unstable manifold are then as follows; here described 
for the case of a Halo orbit in the Hi family. 

1. The libration points, which are the equilibria of Eq. (Q]), are easily determined JSzebehelvll 19671) . In our continuation 
context we note that they have zero velocity components, as well as z = 0, and that for varying /i their x and 
y components lie on connected curves, as shown in Fig. [5] Starting from, for example, the curve of equilibria 
x 2 +y 2 — 1, that exists when ji = (the curve containing the point q in Fig. [5]), the libration points bifurcate from 
x— 1/2, y = ±y/3/2 and y = 0, x = ±1 and we can reach each of the libration points at any given nonzero value of 
/i via a connected path. The eigenvalues of the target libration point(s) are also computed. 

2. Compute the target Halo orbit in the Hi family for which the unstable manifold is to be computed: 

The libration point Jzfi has two pairs of purely imaginary eigenvalues and a pair of real eigenvalues. The two pairs 
of purely imaginary eigenvalues correspond to the planar Lyapunov fa mily Li and V ertical family Vi that bifurcate 
from Jjf\. Compute the family L\, using a standard starting procedure (Doedel ll98ll) at the libration point Jjf\. The 
free problem parameters in this continuation are the period T and the unfolding parameter a. Along Li two branch 
points are located; see Fig. [3] Branch switching at the first of these branch points gives the Halo family Hi. 
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Fig. 5: Computation of the libration points. The point q on the circle of equilibria is a suitable starting point for de- 
termining the libration points for a given nonzero value of jl by continuation. For clarity, the five libration points are 
indicated in the diagram for the case /I = 0.25, rather than for the Earth-Moon case (/i = 0.01215), as /I is close to zero 
in the latter case. The small squares are branch points. 

3. Determine the eigenfunction of a selected Halo orbit: 

Select an appropriate periodic orbit in the Hi family that has one real Floquet multiplier with absolute value greater 
than 1, so that its unstable manifold is two-dimensional. Compute the corresponding unstable eigenfunction as 
follows: Couple the boundary value equations for v in Eq. @ to those for u in Eqs. ©, ©, (|5]). Supplement this 
extended system by an appropriate continuation equation of the form of Eq. ([6]), that also includes v, A, and p, with 
v and p initialized to zero, and A initialized to the desired Floquet exponent as obtained from the decomposition of 
the Jacobian of the collocation equations. Written out, this gives 



In our continuation context we observe that the target period orbit u(f), together with v(t) = and p = (a,A,p) = 
(0, A , 0)), corresponds to a branch point, from which a solution family (u, v; <7, A , p) bifurcates. Along this bifurcat- 
ing family the orbit u(t) remains equal to the target periodic orbit, the Floquet exponent A and the period Tofu(t) 
remain constant, while the unfolding parmameter o remains zero; all to numerical precision. On the other hand, the 
Floquet eigenfunction \{t) becomes nonzero, and consequently p also becomes nonzero. In fact, one continuation 
step along the bifurcating family would be enough to obtain the nonzero eigenfunction v(f ). However, for numerical 
purposes we typically do a few continuation steps along the bifurcating family until p is equal to 1, i.e., until the 
norm of the Floquet eigenfunction equals 1 . This simple procedure to determine the Floquet eigenfunction fits well 
into the continuation and branch-switching algorithms in AUTO, it is numerically stable, and has the additional ad- 
vantage that the Floquet eigenfunction, once computed, can be continued with fixed nonzero p and varying energy. 



u(0 = rf(u(0,cr), 

11(1) =11(0), 




(9) 



v(l)=v(0), 
(v(0),v(0))=p, 
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In that case p = (a, A, T) in the last constraint of Eq. ([5]). This allows the determination of the eigenfunction v in 
highly sensitive cases, for example, for very large or very small Floquet multipliers. 

Once the periodic orbit and its appropriate eigenfunction have been computed, we can proceed to compute the manifold. 
4 Basic Computation of Unstable Manifolds 

Here we describe a method for computing unstable manifolds based on the continuation of orbits that lie in it. These 
orbits start (as a function of time) in the linear approximation of the unstable manifold that was computed above, and 
end in a section where one of the coordinates is fixed. Other possibilities to constrain the end point include fixing the 
integration time or the arclength; we found the fixed end section to b e the most appropriate here. A H these approaches 
are robust against sensitive dependence on initial conditions; see also lKrauskopf and Osingal J2007I) . 




Fig. 6: Plot of a periodic orbit u(f) having Floquet eigenfunction v(f). Also shown is the orbit r(f) in the unstable 
manifold that starts at the point r(0) = u(0) + £iv(0) and ends in a section E. Here u(t) is an Hi Halo orbit, where 
T — 2.3200, E — —1.5052, A = 1.4534, e\ — 0.05, and the section E is at x — 0. For accuracy the absolute value of £\ 
should be smaller, but for clarity of the diagram it is taken larger here. 

We now explain the procedure for computing the unstable manifolds of a given periodic orbit u(f ) in some detail. 
We assume that u(t) has one real Floquet multiplier with absolute value greater than 1, with associated eigenfunction 
v(f), computed as described in Sect. [3] The starting data needed are a point on the periodic orbit, namely u(0), and the 
corresponding value of the Floquet eigenfunction, namely v(0). For a given appropriately small value of £i, take the 
point r(0) = u(0) + £iv(0) as the starting point of an orbit r(t) in the unstable manifold; see Fig. [6] Here t denotes the 
time along the orbit r(f). Similar to the case of the periodic orbit u(t) and its eigenfunction v(f), where for numerical 
reasons the time interval was rescaled from [0, T] to [0, 1], for the non-periodic orbit r(f ) we also rescale time to the unit 
interval. Select a section E, for example, at xz =0 or xz = —0.25, where the orbit r(t) is to terminate; i.e., r(l) £ E. 
(We may allow the orbit r(t) to cross E several times, before it actually terminates in E.) The part of the manifold to be 
approximated is then given by the set of orbits 

{r(0 |r(0)=u(0) + £v(0) and r(l)eZ, for e x < e < e 2 } . (10) 

The range of values of £, namely, [£i,£2), should be chosen to correspond to a fundamental domain, see Fig [6] This 
ensures that the full manifold is swept out, at least locally near the periodic orbit, as £ is allowed to vary from e\ to £2. 
If the interval [£1, £2) is chosen too small then the manifold would be incomplete. Taking this interval too large would 
lead to duplication of orbits, albeit having different lengths. A fundamental domain is such that the orbit that starts at 
u(0) + £iv(0) closely passes the line given by u(0) + £v(0) again, for the first time, at u(0) + £2v(0). For a given value 
of £1, using fundamental properties of the eigenfunction v(f) we can compute the corresponding linear approximation 
of £2 using £2 = e x £\. 

We also note that if £1 is too small in absolute value then the orbit r(t) remains close to the periodic orbit for a long 
time before it escapes to ultimately reach the section E. If, on the other hand, £1 is too large in absolute value then the 
linear approximation of the manifold is no longer accurate. 
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The boundary value problem for r(f ) is then given by 



r(0 = r r f(r(0,0), 
r(0) = u(0)+£v(0) , 

r(l), = x L , (ID 

/ (r(0 - r (t)J (t))dt + (e - £ )£o + (T r - T r0 )T r f = As , 
Jo 

where o = in f(r(f), ct), and E denotes the section x = X£. As done earlier for the periodic orbit and for its eigen- 
f unction, the actual integration time T r of the orbit r(t) appears explicitly in the differential equation, due to the above- 
mentioned scaling of the time t. The pseudo-arclength constraint, the last equation in Eq. dTTb plays a crucial role in the 
algorithm. For a given step size As the pseudo-arclength constraint ensures in particular that the next computed orbit 
r(t) is close to ro(t) over the whole trajectory, and that also T r is close to T r0 . The step size As could in principle be 
constant, but for efficiency and robustness it is adjusted after each successful continuation step, depending on the speed 
of convergence of the Newton/Chord iterations. 

Given the point u(0) on the periodic orbit, and the point r(0) = u(0) + £iv(0) in the direction of the unstable 
manifold at time t = 0, the complete procedure to compute the manifold is then as follows: 

1. Compute a starting orbit in the manifold, using continuation to do the "time integration". More precisely, we use 
numerical continuation with T r as free parameter to compute a "family" of solutions, i.e., the same trajectory, but for 
a set of increasing integration times T r . The equations used are 

r(0 = r r f(r(t),0), 
r(0) = u(0)+£iv(0), 

which correspond to Eq. dTTb , but without the end point constraint. Note that £ is fixed at £ = £\ in each continuation 
step of this starting procedure. The continuation is stopped when r(l) intersects the plane E. (This termination point 
need not necessarily be the first such intersection.) The starting orbit in this continuation is the constant solution 
r(t) = u(0) + £iv(0), with T r — 0. Although it may appear inefficient to do a time-integration by continuation, this 
approach has the advantage that it fits very well into the continuation framework of the algorithms in this paper. 

2. Given an orbit that ends in E, as computed above, the unstable manifold is then approximated by further continuation 
of this orbit, now using the boundary value problem in Eq. dTTb , with the end point r(l) constrained to remain 
in E, and with £ and the integration time T r allowed to vary. If £ varies along a full fundamental domain, the 
continuation would sweep out the full unstable manifold, limited only by the termination condition. However, the 
pseudo-arclength constraint in Eq. dTTb limits the size of the change in the orbit and in the the parameters in any 
continuation step. Hence the full unstable manifold is only obtained if no "obstacles" are encountered, where for 
example T r goes to infinity. As will be seen in the next section, these obstacles can in fact be of much interest. 



5 Example Computations of Unstable Manifolds and Connecting Orbits to Invariant Tori 

In this section we illustrate the basic boundary value technique described above by computing unstable manifolds of 
the Vertical family Vi and the Halo family Hi. For certain ranges of their energy (the thin solid curves in Fig. [2]) and 
period these families have periodic orbits with one real, positive Floquet multiplier outside the unit circle, so that their 
unstable manifold is indeed two-dimensional. The algorithm applies in principle also to higher-dimensional unstable 
manifolds, by using a selected unstable Floquet multiplier, typically the largest one, and its associated eigenf unction. As 
already mentioned, the algorithm also applies to stable manifolds. However, here we restrict attention to examples with 
two-dimensional unstable manifolds. 

Fig. [7] shows the computed unstable manifold of an orbit from the Vertical family Vi . Here the section E is taken at 
xz = 0. As seen in the figure, the orbits in the manifold terminate in E at their second intersection with this plane. Note 
also that the intersection curve of the manifold with the section E has a similar shape as the figure-eight Vertical orbit 
from which the manifold originates. 

Fig. [8] shows the computed unstable manifold of a Halo orbit from the Hi family. The plane E is located at xz — 
—0.25. Note that the manifold changes shape as it propagates. Here, as well as for the unstable manifold of the Vi orbit, 
there is no contradiction in the fact that the cross section of the manifold with the plane E is a self-intersecting curve, 
since we are viewing a projection of the manifold from R 6 (spatial and velocity variables) into R 3 (spatial variables 
only). 

If the section is taken at certain other value-ranges of xz then one may encounter obstacles that prevent the contin- 
uation to cover the full manifold. More specifically, the initial value of the orbits, as determined by the value of £, does 



10 



Fig. 7: The unstable manifold of a periodic orbit from the Vi family. The periodic orbit, labeled Vi, is located on the 
Jzfi side of the Moon, with period T — 3.7700 and energy E = —1.5164. The terminating plane E is located at xz = 0. 
The first intersection of the manifold with E is indicated by the curve labeled E\, and the second intersection, where the 
manifold computation is terminated, is labeled E^. 




Fig. 8: The unstable manifold of a periodic orbit (labeled Hi) from the Hi family. The periodic orbit, which is on the ££\ 
side of the Moon, has period T = 2.5152 and energy E = —1.5085. The terminating plane E is located at xz = —0.25, 
and the first and second intersections of the manifold with E are labeled E\ and 1%. 



11 



Fig. 9: Part of the unstable manifold of the Vi periodic orbit from Fig. [7] together with a superimposed longer orbit in 
this manifold, computed as described in the text. 




Fig. 10: A separate view of the longer orbit from Fig.[9j which winds around a torus near the original Vi periodic orbit. 
The orbit ultimately returns to the plane E that is located at x% =0.5. We remark that in this and subsequent figures 
the orbit color changes from a deep sky blue (RGB code (0,0.5,1)) via gray (code (0.5,0.5,0.5)) to dark orange (code 
(1.0,0.5,0)) as the scaled time increases from to 1. 
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Fig. 11: An orbit obtained by initial value integration for the same value of £, within numerical precision, as used for the 
orbit in Fig. [To] Note that the orbit approaches but then diverts from the torus seen in Fig. [Toj without winding around 
it. 




Fig. 12: Continuation of a longer orbit within the unstable manifold of a Halo orbit with E = —1.5532 with the plane L 
located at xz = 1.02, that is, on the J^-side of the Moon. This orbit winds around a torus near an orbit from the Halo 
family H2, before returning to and ending in the plane E. 
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not cover the entire fundamental domain, but approaches a particular value where r r — )► oo. One possible scenario is that 
the orbit picks up additional loops around a torus-like object near its end point, if such an object exists, before returning 
to and ending in the plane E. This computational phenomenon is not exceptional in the CR3BP; in fact it is easy to find 
specific examples. One such instance is given in Fig. [9] which again shows the unstable manifold of the Vertical orbit 
in Fig. |71 but with a superimposed longer orbit that lies in the same unstable manifold. This longer orbit alone is shown 
separately in Fig.[lOl As is evident from Fig.[lOl the orbit returns to a neighborhood of the original Vertical orbit, where 
it winds around a torus-like object (a quasi-Vertical orbit). 

More specifically, the phenomenon shown in Fig. [9] and Fig. [TCI results from "growing" (see step 1 at the end of 
Sect.HJ a longer initial orbit in the unstable manifold of the Vi periodic orbit, with £ fixed, until it intersects a plane E 
located at xz —0.5. Continuing this orbit with the end point constrained to remain in E, and with £ and T r allowed to 
vary, then results in it winding around a torus near the original Vi periodic orbit, before returning to and terminating in 
the plane E. In contrast, trying to obtain such an orbit directly by shooting techniques appears to be difficult, as such 
orbits approach but then divert from the torus; see Fig.[TT] Indeed, for the orbits to agree it may be necessary to compute 
the BVP orbits to higher precision, i.e., more than the double precision used in AUTO. In particular, this would yield the 
initial value to higher precision. Likewise, the initial value integration, starting from such a more accurate initial value, 
may also need to be in high precision arithmetic. 

A similar result is shown in Fig.[l2lfor an extended orbit from the Hi manifold shown in Fig. [8] Here an orbit in the 
unstable manifold of the Hi periodic orbit is grown until it reaches a plane E located at xz — 1.02, i.e., on the «£?2-side 
of the Moon. Constraining the end point of this longer orbit to remain in E, and allowing £ to vary, then results in the 
end portion of the orbit to wind around a torus near an orbit from the Halo family H2 (a quasi-Halo orbit), as shown 
in Fig. [T2] Meanwhile the initial value r(0) of the orbit approaches a point in the fundamental domain, without fully 
covering that domain. Ultimately the orbit returns to the plane xz = 1.02, as required by the computational set-up. 

Similarly, in Fig.[l3]an extended orbit from the Hi manifold is grown until it reaches a plane E located at xz = 0.6, 
after going around the Earth once. Constraining the end point of this orbit to E, while £ and the integration time T r are 
allowed to vary, results in an orbit that winds around a quasi-Halo orbit near Hi. 

As a final example in this section we compute an initial orbit in the unstable manifold of an Hi Halo orbit, but now 
in the part of the unstable manifold on the side of the Moon. The terminating plane E was taken at xz = 1.02. Further 
continuation of this orbit with xz constrained to remain in E, and with £ allowed to vary, results in the end portion of 
this orbit winding around a quasi-Halo orbit near an H2 Halo orbit, as shown in Fig. [T4j before returning to and ending 
in E. 

In this section we have computed the unstable manifolds of some orbits in the Hi Halo and Vi Vertical families. 
For certain values of £ we find what appear to be heteroclinic connecting orbits from the original periodic orbit to an 
invariant torus. These connecting orbits appear to be more difficult to obtain with shooting techniques. This was seen 
in Fig. [TT] where shooting failed to reveal the invariant torus, when starting near the numerical value of £ for which 
the heteroclinic connection was found in the BVP approach. This is in part due to the fact that the tori are evidently 
unstable, with saddle-type stability since the connecting orbit approaches it, but ultimately also leaves the neighborhood 
of the torus. 

The next objective is to continue the orbit-to-torus connections as the energy E is allowed to vary, and with it the 
base periodic orbit itself. In the following section we explain the computational set-up, while examples are given in 
Sect. |7] In particular we will see in Sect. [7] how such continuation can lead to homoclinic and heteroclinic connecting 
orbits between periodic orbits in the same or different families. 



6 Continuation of Connecting Orbits 

The examples in the preceding section provide evidence for the existence of orbits in the unstable manifold of certain 
periodic orbits that connect to toroidal objects in phase space. These connections exist for specific initial points in the 
fundamental domain, i.e., along the line u(0) +£v(0), £ E [£i,£2). The connecting orbits are generated numerically as 
the initial value of the orbits approaches a final point inside the fundamental domain, as the integration time T r increases 
and additional windings are generated near the end of the orbit. As a consequence the fundamental domain is not fully 
covered, although this can be remedied by repeating the entire sequence of steps in the algorithm, starting from a value 
of £ that corresponds to a point in the part of the fundamental domain that is not covered. 

The question arises as to what happens to the periodic-orbit-to-torus connections if the periodic orbit is varied along 
one of the primary or secondary families. One approach would be to repeat the computations for each one of a sequence 
of periodic orbits along a given family, e.g., along the Halo family Hi. However, here also, continuation is a more 
effective tool that allows interesting connecting orbits to be detected easily along the continuation path. 

We are then led to consider the following approach: Collect the equations that define the periodic orbit (Eqs. ©, 
([5])), the equations defining the Floquet multiplier and eigenfunction (Eq. ©), and the equations for the orbit in the 
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Fig. 13: Continuation of a longer orbit within the unstable manifold of a Halo orbit with E = —1.5631 with the plane E 
located at X£ —0.6. This orbit winds around a torus near the orbit from the Halo family Hi where it started from, before 
returning to and ending in the plane E. 




Fig. 14: An orbit in the Moon-side part of the unstable manifold of an Hi periodic orbit. This orbit is found by fixing 
its end point to remain in a plane E located at x% = 1.02, that is, on the ^-side of the Moon. The orbit winds around a 
torus near an orbit from the Halo family H2, before returning to and ending in the plane E. 



unstable manifold (Eq. (fTTT)). For clarity, this complete set of coupled equations is reproduced in Eq. (fTSt . 
u(0 = rf(u(t),a), 
u(l) = u(0) , 

/ (u(0,u o (0)^ = 0, 
Jo 

v(0 = rf u (u(0,o)v(0 + Av(0, 

v(l) = v(0) , 

\ J \ ) , (13) 
(v(0),v(0)>=p, 
r(0 = r r f(r(0,0), 
r(0) = u(0) + ev(0) , 
r(l)x = x E , 

/ (w(0-w o (0,w , o(0>* + <P-Po,Po> w(0 = (u(0,v(0,r(0), p=(r,a,A,e) 

Jo 



The last constraint shown in Eq. ( fT3t is the suitably expanded version of the pseudo-arclength constraint (Eq. ©) that 
defines the continuation step size. 

Recall that each of the vectors u, v, and r, is six-dimensional. A simple count shows that Eq. ( fT3t represents a system 
of 18 ODEs, subject to a total of 21 boundary and integral constraints, not counting the pseudo-arclength constraint. 
Generically, the continuation of a solution to Eq. (\3[ then requires four free scalar variables. The appropriate choice 
of these parameters is T, a, A, and e, where T is the period of u(f), a is the unfolding parameter, is the Floquet 
multiplier, and £ corresponds to the step size in the direction of the unstable manifold. 

The BVP in Eq. (l3l does not directly define a connection from a periodic orbit to a torus, but a connection from 
a periodic orbit to a section E for a sufficiently large, fixed, value of the integration time T r of r(f). The torus is 
then indirectly continued. This indirect approach is very useful as the d i rect approach is known t o be of considerabl e 
algorithmic compl e xity, as addressed for example in iDieci et al. | dl99lh : iDieci and Loren3 dl995h : lEdoh et alJ dl995h : 
iHenderson I (EoO^ i lSchilder et alJ d2005h : IOlikara I bOlCh . 

The strategy of indirect continuation is somewhat analogous to the computation of a simple homoclinic orbit, i.e., 
an orbit in phase space that approaches a given saddle equilibrium in both positive and negative time. The continuation 
of such a homoclinic orbit in two parameters can be directly formulated as a boundary value problem with asymptotic 
boundary conditions that compute the saddle point and its relevant eigenspaces. In fact, this approach has been used 
very effectively for the continuation of homoclinic and heteroclinic orbi ts, including higher co-d imension orbits, both 
for orbits homoclinic or heteroclinic to equilibria or to periodic orbits (IChampnevs et al. 1996). However, a generic 
homoclinic orbit can also be very effectively approximated indirectly by a periodic orbit of high period, which renders 
the 2-parameter continuation of a homoclinic orbit as simple as the continuation of a periodic orbit (using the techniques 
given in Sect. [3} where the period T is fixed at a sufficiently large value. 

In the next section we demonstrate the use of the indirect periodic-orbit-to-torus approach described above, by 
applying it to the connecting orbits from a Halo orbit in the Hi family to a torus, as initially computed in Sect. [5] We 
will see that these continuation calculations can lead to approximate connecting orbits from the base periodic orbit to 
resonant periodic orbits and to other periodic orbits. The detection of such approximate connections can be refined, 
which we do not do in this paper, to provide more accurate approximations to such connecting orbi ts. These, in turn, 
can th e n be continued directly, usi ng known algorithms for this purpose, as described, for example, in lChampnevs et al] 
(ll996h : lDoedel et al.l d2008L l2009h . and as implemented in AUTO. 



7 Showcasing Three Families of Connecting Orbits 

In this section we describe the results of the further continuation of orbits that connect a periodic orbit to a torus, as 
initially found in Sect. [5] This is done using the 18 -dimensional system in Eq. (TT3T ). which was discussed in the preceding 
section. 

As a first example we start from the orbit shown in Fig.[l3]which connects the Hi Halo orbit with period T — 2.5152 
and energy E — —1.5164 to a quasi-Halo toroidal orbit near the originating Halo orbit, after making one loop around 
the Earth. The continuation procedure using the 18 -dimensional ODE in Eq. (TTTI) allows the energy to change, and with 
it the Hi orbit itself, the connecting orbit, and the torus it approaches. Interesting transitions are encountered along the 
continuation path, namely, in the way the changing connecting orbit winds around the changing torus. 

Examples are shown in Fig. [15] Panel (a) shows the connecting orbit evidently approaching a 5:1 resonant orbit. 
In panel (b) the quasi-Halo orbit has shrunk so it can not be visually distinguished from the periodic Halo orbit in 
whose unstable manifold it lies, that is, the orbit appears to be homoclinic to the Hi orbit. Such an orbit could be 
interesting in space mission design, allowing for occasional large spatial excursions from an Hi orbit for negligible 
energy cost. Panel (c) shows a heteroclinic connection between an Hi orbit and the planar Li orbit with the same energy 
E — —1.5754. Recall that the Li family is the Lyapunov family that bifurcates from the libration point Jzfi, as seen in 
Fig. 12 This heteroclinic connection could be used for a low cost transfer orbit between the Hi and Li orbits with energy 
E — —1.5754, even though these orbits are different in their dynamical properties (see Figs. [3] and© and far apart in 
the bifurcation diagram of the orbits (see Fig. [2}. Panel (d) is similar to panel (b), but now the Halo orbit is connected 
to its corresponding southern Halo orbit, which is the mirror image in the z = of the Hi orbit. That is, the connection 
appears heteroclinic instead of homoclinic. The energies of these connecting orbits and the periods of the periodic orbits 
that they connect to, for these and other connecting orbits in this paper, are given in Table [T] 

During the same continuation one encounters the connecting orbit seen in Fig. [16] There is still ample evidence of an 
underlying torus as this figure shows. However, on this torus the connecting orbit oscillates between a northern Axial Ai 
orbit and its symmetric southern counterpart, each time spending a significant number of rotations very close to each of 
these two periodic orbits. This suggests a low cost transfer orbit between an Hi Halo orbit and the northern or southern 
Ai Axial orbit with the same energy E = —1.5085. The connecting orbit between the northern and southern Axial orbits 
is generic, as will be shown in the next section. Evidence of tori close to heteroclinic connections between two symmetric 
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Fig. 15: A family of connecting orbits from Hi Halo orbits to tori, making one loop around the Earth. The terminating 
plane E is located at xz = 0.6. (a)-(d): Connecting orbits from Hi Halo orbits to (a) a 5:1 resonant orbit near the 
corresponding southern Halo orbit, (b) approximately the originating Hi Halo orbit, (c) an Li planar orbit, and (d) a 
southern Hi orbit, (e): Energy of the Halo orbit versus log(— e) along the continuation path. The continuation terminates 
in the top right corner when the Floquet multipliers of the Hi orbit become complex, and at label 16 in the heteroclinic 
connection shown in Fig. [16] Labels a-d,[l3][T6] and[T9la) and (b) correspond to the values in the panels above and in 
Figs. [T3] [T6j and [19] (f): The logarithm of the unstable Floquet multiplier of the Halo orbit versus log(— e) along the 
continuation path. 
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(a) 




(b) 



Fig. 16: (a): The unstable manifold of the Hi periodic orbit from Fig. [8) together with a superimposed longer orbit that 
was found by continuation using the 18 -dimensional system in Eq. (fT3l) . (b): A separate view of the orbit in the unstable 
manifold. This orbit returns near the originating Hi periodic orbit, where it winds around a quasi-Halo orbit, while 
spending significant time near a northern Axial Ai orbit and its symmetric southern counterpart. As is clearly visible, 
and as required by the computational formulation, the orbit ultimately returns to the plane L, located at xz — 0.6. 
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Axial families can also be found in Figure 1 of lGomez and Mondelol ( 2001 ) for an energy value (corresponding to Fig. [19 



discussed later) not far from that of Fig. [TSJ In Gomez and Mondelol J200 ll) the Axial orbits are denoted by diamond- 
shaped fixed points. 

Fig. [T5le) depicts the change of £ versus the energy, where the labels correspond to the values of panels (a)-(d) and 
Figs.[T3j|T6j and[l9] Note that for most energy values in this range there exist two connections with different values of e. 
The continuation terminates at two points on the right hand side of this diagram, because the Newton-Chord method that 
AUTO employs no longer converges there. One of these points corresponds to the special connecting orbit in Fig. [16] For 
the other termination point the unstable real Floquet multiplier reaches the unit circle, so the two-dimensional unstable 
manifold ceases to exist. Fig.[T5[f) and the thin-to-thick curve transition for Hi in Fig.[2]a) show this Floquet multiplier 
behavior. 

The apparent connections from an Hi Halo orbit to an L2 planar orbit, an H2 Halo orbit, and 5:1 and 6:1 resonant 
orbits shown in Fig. fT7l a)-(d) result from the continuation of the connecting orbit in Fig. [12] again using the 18- 
dimensional system in Eq. dTTb . In this case the connecting orbit makes four loops (instead of just one as in Fig. [15]) 
around the Earth. Here also, the continuation passes through many resonances. Specifically, the connecting orbits shown 
in panels (c) and (d) of Fig. fT71 approach a period-5 orbit and a period-6 orbit, respectively. 

We reiterate that the originating Hi Halo orbit changes during the continuation, and with it its period and energy. 
Consequently the energy of the connecting orbit as well as the energy of the torus that it approaches change with it, as 
these energies are identical. Similar to panels (e) and (f) of Fig. [15] Fig. [TTl show diagrams for the energy and Floquet 
exponents. Note however, that in this case the continuation does not terminate but loops around, and for every energy 
value between the two extrema there exist two connecting orbits. Numerically it was observed that the value of log(— e) 
is always close to — 2A, so to visualize the loop these two quantities were subtracted in panel (f). 

Note that the two extrema in panels (e) and (f) are detected as folds by the continuation software, and that these 
could in turn be followed, for instance by adding the mass-ratio /1 as a free parameter. The folds themselves do not 
appear to correspond to particularly interesting orbits (panels a and b in Fig.[T7]are close to but not at the fold locations). 

Similarly, Fig. [T8l a) shows an interesting orbit from the continuation of the connecting orbit in Fig. [14] namely a 
direct connection, without looping around the Earth, from an Hi orbit to an L2 Lyapunov orbit. Fig.fTOb) and (c) show 
diagrams much like those in Fig. [T7] where the continuation curve is a loop. However, compared to Fig. [T7] here the 
range of the energy level and the differences between values of log £ are greater. Also note that the sign of £ is always 
positive, rather than negative as before, as we are now considering the Moon side of the unstable manifold of the Hi 
orbits and so the sign of £ in Eq. (\3[ is changed. 



8 Existence of Connecting Orbits 

As seen in the preceding sections, boundary value formulations provide a powerful tool to compute unstable manifolds 
and connecting orbits in the CR3BP. Our exploration shows only the tip of the iceberg of a wealth of interesting orbits. 
This section discusses how the connecting orbits that were found relate to the existing literature. As mentioned in Sect.Q] 
the periodic orbits and tori themselves, as well as homoclinic and heteroclinic orbits connecting Li and L2 periodic orbits 
have been studied extensively. However, for the spatial case, to the best of our knowledge, connecting orbits from Hi 
Halo to quasi-Hi or quasi-H2 hav e not been explicitly found before. We must mention that connections from quasi-Hi 
to quasi-H2 orbits can be found in lGomez et all J2004I) . 

In all three continuations in Sect. [7] a collection of interesting objects was seen, including apparent heteroclinic 
connections between Halo orbits and resonant orbits, i.e., seemingly heteroclinic connections between Halo orbits and 
/z-periodic orbits, where the complex Floquet multipliers of the Halo orbit near the ^-periodic orbit are close to e 1%l l n . 
Numerical data corresponding to figures in this paper are given to 5 significant digits in Table [T] 

The existence of certain connecting orbits can be explained by counting dimensions. For Fig. [16] the multipliers of 
the symmetric Axial orbits give rise to three-dimensional stable and unstable manifolds, so the connection between the 
two northern and southern Ai Axial orbits seen in Figure [T6l is generic for the CR3BP posed in R 6 . The connection 
from the Hi Halo orbit to the Ai Axial orbit is codimension-one, since the Hi Halo orbit has a two-dimensional unstable 
manifold, and the stable manifold of the Ai Axial orbit is three-dimensional. Similarly the planar Li and L2 Lyapunov 
orbits for the energy values that we observe in Figs.[l5]c),[T7]a), and[THa) have three-dimensional stable and unstable 
manifolds, so connections from Hi Halo Orbits to those orbits are also codimension-one. By standard bifurcation the- 
ory, codimension-one connecting orbits can be expected to arise as one parameter is varied, for specific values of that 
parameter. This is indeed what we observed in these cases, where we found such connecting orbits for specific values of 
the energy E. 
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Fig. 17: A family of connecting orbits from Hi Halo orbits to tori, making four loops around the Earth. The terminating 
plane L is located at xz — 1.02. (a)-(d): Connecting orbits from Hi Halo orbits to (a) a planar L2 Lyapunov orbit, (b) 
approximately an H2 Halo orbit, (c) a 5:1 resonant orbit near the libration point j£?2, and (d) a 6:1 resonant orbit, (e): 
Energy versus log(— e) along the continuation path. The continuation curve is a loop. Labels a-d and fT2l correspond to 
the values in the panels above and in Fig. [12] (f): The logarithm A of the unstable Floquet multiplier of the Halo orbit 
versus log(— e) + 2A along the continuation path. Here log(— e) + 2A is plotted instead of log(— e) to make the loop 
visible. 
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Fig. 18: (a): An orbit directly connecting the Moon side of the unstable manifold of the Hi orbit with E — —1.5728 
to an L2 Lyapunov orbit, obtained by continuation of the orbit shown in Fig.[l4] (b): Energy versus log(e) along the 
continuation path. The continuation curve is a loop. Labels a and fT4l correspond to the values in the panels above and 
in Fig. [14] (c): The logarithm of the unstable Floquet multiplier of the Halo orbit versus log(e) along the continuation 
path. 



8.1 Connecting Orbits from a Halo orbit to a quasi-Halo orbit 



In contrast to the connecting orbits considered above, the connections from the Hi Halo orbit to the same or other Halo 
orbits cannot simply be explained by counting dimensions. The Hi Halo orbit has a two-dimensional unstable manifold 
and a two-dimensional stable manifold, so such connections would be codimension-two and would not be expected 
along a family (in the sense of this paper) of Halo orbits. The apparent Halo orbits in the Halo to Halo connecting orbits 
may in fact be very thin quasi-Halo orbits (tori). However, even if this distinction were mathematically important, it 
would resumably be of little practical importance in space mission design. 

We justify the presence of the connecting orbits from Hi Halo orbits to quas i -Halo orbits by reca lling the detailed 
analys is of center manifolds around collinear libration points in IJorba etal.Ul999h J Gomez et al. and lKoon et al.l 

J2008I) . For a prescribed energy level, the thin solid curves in Fig. [2 show that there exists a Halo orbit near J^-, for 
i — 1,2, that has two complex Floquet multipliers on the unit circle and different from one. There are two more real 
Floquet multipliers, one of them greater than one and the other less than one. The Floquet multipliers and corresponding 
energy levels for the Halo orbits in Figs. [8] and [12] to [19] are shown in TableQ] 
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Figure 


Energy 


Type 


Period 


Non-trivial Floquet exponents 


I7I9I10I11I 


-1.5164 


Vi 


3.7700 


±6.4948, ±0.077175 -2ni 


m 


-1.5532 


Hi 

quasi-H2 


2.7873 
3.2753 


±6.1730, ±0.21859 -2tt/ 
±5.9348, ±0.17737 -2tt/ 


El 


-1.5631 


Hi 


2.7821 


±6.6179, ±0.16367 -2ni 


m 


-1.5733 


Hi 

quasi-H2 


2.7716 
3.3805 


±7.0356, ±0.11217 -2tt/ 
±6.7839, ±0.068079 • 2ni 


Eta) 


-1.5552 


quasi-Hi (5:1 res) 


2.7868 


±6.2673, ±0.20690 -2tt/ 


Etb) 


-1.5716 


Hi 


2.7736 


±6.9682, ±0.12050 -2tt/ 


Etc) 


-1.5754 


Hi 
U 


2.7690 
2.8982 


±7. 1183, ±0.10189-271/ 
±0.29251, ±7.4268 


Etd) 


-1.5617 


Hi 


2.7832 


±6.5556 ,±0.17135 • 2ni 


I5ITSI 


-1.5085 


Hi 

Ai 


2.5152 
4.0117 


±2.8541, ±0.38928 -2%i 
±0.43035, ±6.1278 


Eta) 


-1.5276 


Hi 

L 2 


2.7472 
3.9550 


±4.7666, ±0.39460 -2%i 
±5.9229, ±0.52353 


Etb) 


-1.5679 


Hi 

H 2 


2.7776 
3.3567 


±6.8207, ±0.13870-2;n 
±6.5802, ±0.095719-271/ 


Etc) 


-1.5349 


Hi 

quasi-H2 (5:1 res) 


2.7715 
3.1170 


±5.2225, ±0.33730-2tt/ 
±4.9006, ±0.30348 • 2ni 


Etd) 


-1.5303 


Hi 

quasi-H2 (6:1 res) 


2.7581 
3.0568 


±4.9440, ±0.23394-2tt/ 
±4.5639, ±0.34345 • 2ni 


Eta) 


-1.5728 


Hi 

quasi-H2 


2.7722 
3.3784 


±7.0158,±0.11462-2tt/ 
±6.7651, ±0.070727-271/ 


m 


-1.5130 


Hi 


2.6176 


±3.5327,±0.45736-2tt/ 



Table 1: Numerically computed values for all figures. The first line in each box contains the data for the starting periodic 
orbit u(^). If the orbit r(^) connects the periodic orbit to a different periodic orbit, or to a torus that surrounds a different 
periodic orbit (a quasi-periodic orbit), then the second row contains data for this second periodic orbit. There are always 
two trivial Floquet exponents equal to zero. 



The unstable manifold of such a Halo orbit is a two-dimensional surface, see Fig. [8] On the other hand, close to 
either of the two libration points for j — 1,2, there is a four-dimensional center manifold, which exists due to the 
fact that Jzf/ has four eigenvalues on the unit circle and two more on the real line excluding one. 

On the energy surface the center manifold of the co-linear libration points is a little different. In lKoon et al.1 J2008b 
the authors verify that for fixed energy, the center manifold of j£fy, for j = 1,2, is a normally hyperbolic invariant 
manifold (NHIM) corresponding to a normally hyperbolic three-sphere that is invariant for the linearized system. Thus, 
in a small neighborhood of Jzf/ the center manifold becomes a deformed three-sphere for the nonlinear system. Moreover, 
it is a three-dimensional hyperbolic invariant manifold in a five-dimensional energy surface with one stable direction. 
Therefore the stable manifold of the center manifold is four-dimensional. 

In this setting, dimension counting gives a heuristic argument that motivates why connecting orbits from Halo orbits 
to torus-like objects appear quite naturally in our problem, as seen in Figs.[l2l[T3l[l4l and [19] We have mentioned that 
the center manifold of a libration point Jz^ restricted to the energy surface is a three-dimensional NHIM, so we can do a 
dimension counting analysis similar to those above. The main observation is th at for a fixed value o f the energy, these 
quasi-Halo orbits are lower dimensional whiskered quasi-periodic solutions (see lFontich et al.l J2009h ) that lie inside the 
three-dimensional center manifolds of the collinear libration points. 

Now, the Halo orbit is a one-dimensional normally hyperbolic object with one unstable direction. The Halo orbit 
itself also belongs to the center manifold of Jz^ and its unstable manifold is a two-dimensional object in the energy 
surface. Finally, by dimension counting, we notice that a transversal intersection between the unstable manifold of the 
Halo orbit of Jzf/ and the stable manifold of the center mani fold of Jzfy is a one-dimensional object. 

Several authors dJorba etal.| [l999: Masdemont 2005; Ale ssi et al.1 12009) have explored the center manifolds of 
collinear libration points by means of Lindstedt-Poincare expansions of the Hamiltonian restricted to each center man- 
ifold. These computations suggest that for the energy levels under consideration in this paper there exist large regions 
in the center manifold that exhibit quasi-periodic motions. The z — sections of the center manifolds reveal that most 
of these quasi-per iodic obits are two-dimensional invariant tori in the center manifold (see for instance Figure 3 in 
I Ales si et all J2009h ). These invariant tori are also normally hyperbolic and their corresponding stable manifolds are two- 
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Fig. 19: (a) and (b): Two orbits for the family displayed in Fig. 1131 where E — —1.5130. (c) and (d): Their intersec- 
ti ons with the plane z — 0. The energy value corresponds to Figure 1 in G omez and Mondelol J200lh and Figure 3 
in lAlessi et al.1 d2009h . where the energy is denoted as h — E+ ja(1 — fi)/2 - 
intersections of the Ai Axial orbits, as in Gomez and Mondelo (200l|). 



-1.507. The two diamonds correspond to 



dimensional objects inside the four-dimensional stable manifold of the center manifold. These codimension-one stable 
manifolds separate the four-dimensional stable manifold into regions. 

A trajectory in the intersection of the unstable manifold of the Halo orbit and the stable manifold of the center 
manifold has a positive probability of being either on the stable manifold of an invariant torus or in a region bounded by 
stable manifolds of invariant tori. It is even possible that the trajectory belongs to the stable manifold of a periodic orbit 
(e.g., a Halo orbit) in the center manifold of J*f). However, this is unlikely since the stable manifold of a periodic orbit 
is a two-dimensional object inside a four-dimensional stable manifold. Thus, it is more likely that a trajectory in the 
intersection of the unstable manifold of a Halo orbit with the stable manifold of the three-dimensional center manifold 
of a libration point passes close to a two-dimensional hyperbolic quasi-periodic orbit in the center manifold. 

Fig. [I9]a) and (b) show two connections from a Halo orbit to a quasi-Halo orbit nearby, from the continua- 
tion displayed in Fig . [T5l where the energy value E — —1.513 matches the value h — —1.507 used for Figure 1 in 
iGomez and Mondelol (1200 ll) and Figure 3 in lAlessi et all (120091) . where h = E + ^(1 + fi)/2. In Fig. QHc) and (d) we 
show t he intersections of these traje ctori es with the plane z = 0. If we compare these intersections to the corresponding 
ones in lGomez and Mondelol (l200ll) and lAlessi et al.l J2009I) . we discover that the trajectory seems to approach one of 
these quasi-Halo orbits and then drift away from it. We conclude that the trajectory computed by our methods is shad- 
owing a connecting orbit in the intersection of the unstable manifold of the Halo orbit and the stable manifold of the 
corresponding normally hyperbolic lower dimensional torus. 
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9 Numerical Aspects 



In our computations we used the continuation and bifurcation software AUTO fooedellll98ll : boedel et al.lEoTob for 
computing families of periodic orbits, associated Floquet eigenfunctions, unstable manifolds, and connecting orbits. 
Orbits are continued in AUTO as solutions to a suitable boundary value problem (BVP), as described in this paper. To 
compute approximate soluti ons of BVPs, AUTO uses the method of Gaus s -Legendre collocation with piecewise poly- 
nomials on adaptive meshes dde Boor and Swartzll 19731 : lAscher et al.|[l98ll) . These calculations can be fast; for example 
the entire Vertical family of periodic orbits Vi can be computed in 0.12 seconds on a dual core laptop, using as few 
as 20 mesh intervals with 4 Gauss collocation points each, and using 48 continuation steps. In our actual computations 
we use more mesh intervals to ensure high accuracy, e.g., with 100 mesh intervals the family Vi can be computed in 
0.25 seconds. We also often use more continuation steps, for example for the purpose of generating data for computer 
animations. 

Two-dimensional unstable manifolds computed as a solution family by continuation, as done in Sect. [4] require more 
mesh intervals and continuation steps when the orbits in the family wind many times around a torus. For example, a 
connection from the Halo orbit Hi to a torus near H2, winding around the torus 12 times, can be reached easily in 15 
seconds, using 200 mesh intervals and 1000 continuation steps. However, we also computed such connections with a 
much higher number of windings, which requires a correspondingly higher number of mesh intervals and continuation 
steps. In fact, we have used up to 1500 mesh intervals in such cases. Furthermore, to ascertain the correctness of our 
results, we have computed these manifolds with various choices of the number of mesh intervals. 

Similarly, the continuation of solutions to the 18 -dimensional system to follow periodic-orbit-to-torus connections, 
as given in Sect. [6] and used in Sect.[7J requires a correspondingly high number of mesh intervals. Since the dimension 
of the system is then 18, i.e., three times the dimension of the systems used for continuing periodic orbits and unstable 
manifolds, the computation time would increase by a factor 3 3 , that is 27. However, this is reduced since AUTO takes into 
account the zero structure of the submatrices of the full Jacobian matrix that correspond to individual mesh intervals. 
On the other hand, the connecting orbit requires significantly more mesh intervals then the base periodic orbit and 
its Floquet eigenfunction, and AUTO does not take advantage of this. Thus the continuation of the coupled system 
(periodic orbit, Floquet eigenfunction, connecting orbit) in Sect.[6]can take significant computer time, also because such 
continuation with varying energy requires many continuation steps. In extreme cases the calculations have taken up to 6 
hours computer time. 
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